Neutron matter at zero temperature with auxiliary field diffusion Monte Carlo 
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The recently developed auxiliary field diffusion Monte Carlo method is applied to compute the 
equation of state and the compressibility of neutron matter. By combining diffusion Monte Carlo for 
the spatial degrees of freedom and auxiliary field Monte Carlo to separate the spin-isospin operators, 
quantum Monte Carlo can be used to simulate the ground state of many nucleon systems (A < 100). 
We use a path constraint to control the fermion sign problem. We have made simulations for realistic 
interactions, which include tensor and spin-orbit two-body potentials as well as three-nucleon forces. 
The Argonne v' 8 and v' 6 two nucleon potentials plus the Urbana or Illinois three-nucleon potentials 
have been used in our calculations. We compare with fermion hypernetted chain results. We report 
results of a Periodic Box-FHNC calculation, which is also used to estimate the finite size corrections 
to our quantum Monte Carlo simulations. Our AFDMC results for v§ models of pure neutron matter 
are in reasonably good agreement with equivalent Correlated Basis Function (CBF) calculations, 
providing energies per particle which are slightly lower than the CBF ones. However, the inclusion 
of the spin-orbit force leads to quite different results particularly at relatively high densities. The 
resulting equation of state from AFDMC calculations is harder than the one from previous Fermi 
hypernetted chain studies commonly used to determine the neutron star structure. 
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I. INTRODUCTION 

The important role played by N-N correlations on sev- 
eral properties of dense and cold hadronic matter is a 
well established fact Q|. Less established are quantita- 
tive studies performed with realistic nuclear interactions 
derived from N-N data and the spectra of light nuclei. 
The strong repulsion at short range accompanied with 
the strong spin-isospin dependence, make ab initio calcu- 
lations of the nuclear matter equation of state one of the 
most challenging problems in strongly correlated many- 
body theory. 

A theoretical calculation of the nuclear matter en- 
ergy per particle, as a function of the number density 
p, the temperature T and the neutron-proton asymme- 
try a = (N — Z)/(N + Z), with an uncertainty of less 
than an MeV has become a fundamental issue. On one 
hand, one would like to use the observational data from 
neutron stars and supernovae, as well as from heavy-ion 
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collisions to get information on the many-body nature of 
the nucleon interaction. On the other hand, it is of in- 
terest to understand the effect of N-N correlations, and 
particularly of those induced by the tensor force, on the 
structure and the evolution of compact astrophysical ob- 
jects iiSili. 

In this paper we limit ourselves to non-relativistic 
model hamiltonians. Modern two-body potentials 0, IE 
9] fit the Nijmegen N-N data [13 below 350 MeV at a 
confidence level of x 2 /^data ~ 1, and to a large extent 
give equivalent results for several nuclear and neutron 
matter properties However, it has become evident 
that a two-body potential alone is not sufficient to re- 
produce the experimental data of nuclei other than the 
deuteron (A — 2). In the last few years, the Urbana- 
Argonne collaboration has produced three-body force 
models which, when added to the two-body potential, 
provide a satisfactory fit to the binding energies and the 
low-lying states of light nuclei with A < 10 [3 [HI CI • 

It would be desirable to have microscopic calculations 
of the equation of state of nuclear matter with an accu- 
racy comparable to that of light nuclei or, at least, on the 
order of the experimental uncertainties of the equilibrium 
density, po, binding energy per particle at po and com- 
pressibility. This can be considered as the minimal re- 
quirement to attempt the study hadronic matter at den- 
sities larger than po, and/or with large asymmetries (a 
close to 1) in a realistic way. Such calculations have to 
deal necessarily with potentials which are strongly spin- 



2 



isospin dependent and which include a three-body force. 

Most of the microscopic calculations of the nuclear 
matter equation of state carried out in the last decades 
have been performed by using perturbation theories 
based either on ladder diagram summation, like Brueck- 
ner or Green's Function theories 0, , or Correlated 
Basis Function theories, based on Fermi hypernetted 
chain techniques 0, 0, . In spite of the important 
advances made in recent years in the above theories, the 
required accuracy for the equation of state has not yet 
been reached. 

Quantum Monte Carlo methods have been very suc- 
cessful in calculating the properties of strongly interact- 
ing systems in condensed matter physics. They are sub- 
stantially exact, apart from statistical errors, finite size 
effects and the well known sign problem [l^ for Fermi 
systems. They have been recently used to p erform quan- 
tum simulations of light nuclei [3, I20I l2l| with modern 
nonrelativistic Hamiltonians of the type discussed above. 
However, the exponential growth in the number of spin- 
isospin states with the number of nucleons A, has kept 
this method from being applied to larger nuclear systems. 

Auxiliary field diffusion Monte Carlo H3 (AFDMC) 
has been especially developed to tackle the problem of 
computing the binding energy of a relatively large nu- 
clear system at the required accuracy. In this approach 
the particle coordinates are propagated as in standard 
diffusion Monte Carlo. Auxiliary fields are introduced 
to uncouple the spin-dependent interaction between par- 
ticles by means of a Hubbard-Stratonovich transforma- 
tion. The particle spins only interact with the auxiliary 
fields which, when integrated, produce the original in- 
teraction. The method consists of calculating the auxil- 
iary field integrations by Monte Carlo sampling and then 
propagating the spin variables. This propagation results 
in a rotation of each particle's spinor governed by the 
sampled values of the auxiliary variables. The result is a 
sampling of the spin variables which should have less vari- 
ance than a direct approach where the spins are flipped. 

The tensor force couples the spin configurations with 
the orbital angular momentum so that the wave function 
becomes complex. The resulting fermion phase prob- 
lem is handled by applying a path-constraint approxi- 
mation analogous to the fixed-node approximation. The 
AFDMC method for the spin-isospin calculations can be 
viewed as a generalization of the method of Zhang et 
al. [23L |24| used in condensed matter lattice systems to 
the spin-isospin states of nucleon systems, while retain- 
ing standard diffusion Monte Carlo for the spatial degrees 
of freedom. The AFDMC method has proved to be effi- 
cient in dealing with large nucleon systems interacting via 
semi-realistic potentials [22, H^, |2(| and spin-polarized 
neutron systems [2?} ■ 

The aim of this paper is to give a detailed descrip- 
tion of the AFDMC method and to report results for 
the equation of state of pure neutron matter (a = 1) 
with a fully realistic nuclear interaction, at zero temper- 
ature. It presents results of AFDMC simulations of 14, 



38, 66 and 114 neutrons in a periodic box, interacting 
via a realistic potential which includes two-body tensor 
and spin-orbit components, as well as three-body forces. 
Particular attention is paid to the 14 neutron system, 
which may serve as a homework problem for different 
many-body techniques. It is small enough to be han- 
dled by traditional quantum Monte Carlo methods |28| . 
However, it will be shown that the finite size effects of 
14 neutron systems are hard to estimate in a realistic 
way. Actually results obtained with larger systems (66 or 
114 neutrons) show that the equation of state of neutron 
matter cannot be simulated starting from 14 neutron in 
a box, particularly in the high density region. Finite size 
effects for the larger systems considered here can be fairly 
well estimated by the rec ently developed Periodic Box 
FHNC (PBFHNC) theory [SjjJ. We have also performed 
AFDMC calculations of the binding energy of symmetric 
and asymmetric nuclear matter. A few results obtained 
with semi-realistic spin-dependent central potentials are 
presented and discussed. 

The plan of the paper is the following. The Hamil- 
tonian used in this work is shown in the next section. 
In section IIHI the problem of the spin degrees of free- 
dom in quantum monte carlo simulations is discussed. 
Section \N\ is devoted to the description of the AFDMC 
method, including the calculation of the spin-orbit and 
the three-body terms of the Hamiltonian. A discussion of 
the finite size effects along with the Periodic box FHNC 
method is given in Section [3 The results for the neu- 
tron matter equation of state are presented and discussed 
in Section IVII The conclusions and perspectives for the 
present work are in Section IVIII 

II. THE HAMILTONIAN 

We use a non relativistic Hamiltonian of the form 

H = T + V 2 + V 3 

= E V " + E^ + E (!) 

] = 1.N j<k j<k<l 

containing the kinetic term, where we have used 
h 2 /{2m) = 20.73554 MeV fm 2 (which corresponds to the 
n—p reduced mass), and two- and three-body potentials. 
The two-body potential belongs to the Urbana-Argonnc 
vi type 

£ 

n = E "i* = E E v P (rjk)0^ ( i; k) , (2) 

j<k j<kp—l 

where j and k label the two nucleons, Tj k is the distance 
separating the two nucleons, and the spin-isospin depen- 
dent operators O p (i, j) for p = 1, 8 are given by 

O p=1 ' 8 (j, k) = (l, a, ■ <J k ,S ]k ,L jk ■ S jk ) <8> (1, tj ■ T k ) ,(3) 
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where Sjk — i{fjk'^j){fjk'^k) — Sj-Sk is the two-nucleon 
tensor operator, and Lj k and Sjk are the relative angular 
momentum and the total spin, given by 

Ljk = J^j-h) x (V,-- VO , (4) 

Sjk = |(<?j +&k) • (5) 

The full Argonne V\% potential consists of €=18 com- 
ponents. Besides the 8 components given in Eq.Q, it 
includes the 6 (L 2 , 1? Bj -a k , (L- S) 2 ) ®(l,Tj -f k ) charge 
independent ones, as well as 4 other charge-symmctry- 
breaking and charge-dependent terms. 

We use a simplified isoscalar version of the vi$ po- 
tential, the so called v' 8 two-body potential [2£|. This 
potential has been obtained with a new fit to the N-N 
data, with only the first eight spin-dependent operators 
in Eq.Q included. It equals the isoscalar part of Vi& in 
all S and P waves as well as in the 3 Di wave and its cou- 
pling to the 3 Si. It has been used in a number of GFMC 
calculations in light nuclei [20f ; as well as FHNC/SOC 
calculations in nuclear matter [l7|]; differences with the 
Vis potential give small contributions and can be safely 
estimated perturbatively or from FHNC/SOC calcula- 
tions. 

For the sake of completeness, we report here the pa- 
rameterization of the Argonne v' s two-body potential. 

8 

v p (r) = J2 A p , m F m (r) , (6) 

m— 1 

where odd and even components refer to r— 
independent and r-dependent operators respectively. 
The spin-independent part vf? of the two-body poten- 
tial is given by the first component v p —i(rij) only. The 
constants A p ^ m and the functions F m {r) are given in the 
appendix. 

In the case of pure neutron matter (PNM), the isospin 
exchange operators are replaced by the identity. 

We denote by v' 6 the two-body potential model ob- 
tained by restricting the v' s potential to the first 6 (3 
for neutron matter) components. Note that this trunca- 
tion of the Argonne v' s should not be confused with the 
recently produced Argonne AV6' potential |3(j . 

The three-body interaction used in our calculations of 
the equation of state is the Urbana IX potential [2(J • For 
neutrons, the Urbana-IX interaction is given by the sum 
of a spin independent and a spin dependent part 

Vjki = VfJ + Vgf , (7) 

where 

yfki = u o ^2 7l2 ( m T)C3;^7)T 2 (m 7r ,c 3 ;r; fc ) , 

cyclic 



V fk? = B 2-k X! i X Jli X lk} ' ( 8 ) 

cyclic 

and the operator XJ k is given by 

x jk = Y(m n ,c 3 ; r jk )S 2 ■ a k + T(m v , c 3 ; r jk )S jk . (9) 

The values of the parameters of the Urbana IX three- 
body potential, used in our calculations are shown in the 
appendix. Notice that in some of our earlier AFDMC 
calculations we have used C3 = 2.0 fm -2 and fj, = 0.7 
fm" 1 , as given in the original papers proposing the Ur- 
bana IX potential [3lJ and the v' a model interaction • 
Changing C3 from 2.0 to 2.1 leads to a ~ 10% additional 
increase of the three-body force contribution in neutron 
matter. In the following, we will denote with AU8' the 
v' 8 plus Urbana IX interaction, with AU6' the v' 6 plus 
Urbana IX interaction. 

We have also considered the recently developed Illinois 
three-body potentials, which include two A intermediate 
state diagrams 13], and denoted with ELI,. . ., IL4. 

III. SPIN DEGREES OF FREEDOM 

Standard Green's function or diffusion Monte Carlo 
methods for central potentials sample only the particle 
positions since the spin or isospin of the particles can be 
fixed. The Green's function Monte Carlo method used 
in light nuclei also samples the particle positions, but a 
complete description of the spin degrees of freedom is 
kept for each position sample leading to an exponential 
growth of the number of spin-isospin states with particle 
number A. This exponential behavior can be removed by 
sampling rather than summing the spin-isospin degrees 
of freedom. 

We define a walker to be the 3A coordinates of the A 
particles and A spinors each giving the four amplitudes 
for a particle to be in the proton up, proton down, neu- 
tron up and neutron down states. For the special case 
where walkers are sampled from the usual neutron-proton 
up-down basis, the spinors would be one of (1,0,0,0), 
(0,1,0,0), (0,0,1,0), and (0,0,0,1) for each particle. 
Our auxiliary field method requires the more general def- 
inition as shown below. 

As usual, the overlap of the walker bra with the trial 
ket is the wave function amplitude, 

(R,S\^ T ) =^t(R,S) . (10) 

Direct sampling of the spin-isospin in the usual spin 
up/down basis requires a good trial function that can be 
evaluated efficiently. This can be most easily seen for the 
variational formalism, but the same analysis applies to 
Green's function or diffusion Monte Carlo. A variational 
Monte Carlo calculation can be formulated by minimizing 
the expectation value of the Hamiltonian, 

(* T \H\* T ) 
[ ' (* T |* T ) 
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J dREs,s> S')H sliS ^ T (R, S) 

where for a vq interaction we would have 



(11) 



Hs>,s = (S'\S) 



2m 



(RS'\V\RS) , (12) 



with a straightforward generalization for spin-orbit 
terms. 

Variational Monte Carlo can be implemented with ei- 
ther spin sums pi Is! Eii] 



(H) = 



P(R) 
E L (R) 



dRE L {R)P(R) , 

J2s\*t(r,s)\z 

fdRJ2s\*T(R,S)\2> 

£ 5S , *t(«. S')H 8 ',s*t(R,S) 



Es\*t(R,S)\< 



(13) 



or spin samples |3£ 



(H) 

P(R,S) 
E L (R,S) 



j dRj2 E L(R, S)P(R, S), 
J s 

\*t(R,S)\ 2 
fdREs\*T(R,SW> 
E s , ** t (R,S')H S ',s*t(R,S) 
\*t{R, S)\ 2 



(14) 



In these equations, P is the probability density to be 
sampled and El is the local energy. A typical variational 
calculation would use the Metropolis algorithm to sample 
either R or R and S from P, and average the value of 
the local energy over these samples. 

Notice that for an eigenstate of H, both El(R, S) and 
El(R) are constant. So, as for central potentials, the 
variance will be low if the trial function is accurate. No- 
tice also that the spin sum S' in the definition of El (R, S) 
is polynomial rather than exponential in A. For example 
a pair potential will have only order A 2 terms where two 
particles have different spin-isospin. 

The variance per sample for complete spin sums will 
be lower than for spin samples. However, since the spin 
sums grow exponentially with particle number spin sam- 
pling should be more efficient for large particle number if 
the trial function can be evaluated efficiently for a single 
many-particle spin state S. 

Unfortunately, all of the good trial wave functions cur- 
rently used for large numbers of particles cannot be eval- 
uated efficiently for a single many-particle spin state S. 
For example light nuclei variational Monte Carlo calcu- 
lations are typically done using a pair product (or more 
complicated) wave function, 



l*P> = 



p 



m, as) 



where S symmetrizes the operator products, and |<&) is 
the antisymmetric model state. While the symmetrizer 
produces all possible orderings of the operators and 
therefore gives 0(A 2 \) terms, normally the commutator 
terms are fairly small and the ordering of the operators 
is sampled. However, even within a fixed ordering, each 
operator in the product term when operating on a sin- 
gle many-particle spin-isospin state will produce 4 or 8 
new states depending on whether isospin exchange gives 
a new state. 0{A) operators out of the 0(A 2 ) total act- 
ing on a single state are enough to populate all the states. 
Therefore a straightforward evaluation of (RS\^fp), for 
this wave function will have the same computational com- 
plexity as evaluating a complete set of spin-isospin states 
at the position R. Since computing all the states have 
the same cost as a single state, full spin sums are used 
for these calculations. 

If good trial functions for spin-isospin dependent inter- 
actions can be devised which can be evaluated or sam- 
pled efficiently at a single many-particle space position 
and spin-isospin state, straightforward generalizations of 
standard central potential Monte Carlo methods, both 
variational and Green's function, with spin-state sam- 
pling will solve the nuclear many-particle Hamiltonian. 



IV. THE AFDMC METHOD 

Since direct evaluation of the pair product wave func- 
tion is not computationally feasible for large numbers of 
particles, and so far we have no good methods of sam- 
pling these wave functions, we instead drop the operator 
terms altogether and sample the spin-isospin variables 
using a rather poor, but easy to evaluate, wave function. 
Since this wave function does not contain amplitudes of 
the spin states of the correct solution we cannot use it 
to sample the spins. Instead, we rewrite the propagator 
as an integral over auxiliary fields using the Hubbard- 
Stratonovich transformation 



\xo At 



1 



dxe 



(16) 



where O can be a one-body operator. To make use of this 
transformation we write our propagator as the left-hand 
side of Eq. QUI so that the integrand of the right hand 
side is a product of one-body terms. The integrand has 
a form such that propagating a walker at \R, S) results 
in another walker of the same form at \R', S'). 

For N neutrons, the v§ two-body interaction can be 
split into two parts 



i'ik 



/ , 3 

j<k 



B 



j,ct,k,{3 



A 



■j,a;k,/3 , 



(17) 



where roman subscripts like j and k are particle labels 
while greek subscripts like a and (3 are cartesian compo- 
nents. The matrix A and the scalar B are functions of 
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the particle positions, 



where 



j<k 

Aj,a;k,/3 = (V3(rjk) + Vi(rjk))S a /3 + 

[vs(rjk) + ve(rjk)] • x a f jk ■ xp - 5 a/3 ] . 

(18) 

Aj.u;k,i3 is taken to be zero when j = k. A can be viewed 
as a 3./V by 3./V real symmetric matrix. It therefore has 
real eigenvalues and eigenvectors defined by 



k,/3 

The potential can be written as 

^Vjk = B+- Yj ^a^n^nlk^. 



(19) 



j<k 



j,a,k,P,n 
3A 



with 



j a 



(20) 



(21) 



Each of the O n is a sum of 1-body operators as required 
above. 

After applying the Hubbard-Stratonovich transforma- 
tion, the short time approximation for the propagator 
can be written as 



1 

W 



3A/2 

exp 

XI 



m\R-R'\ 2 
2h 2 At 



-B(R)At 



(22) 



The O n do not commute, so we need to keep the time 
steps small so that the commutator terms can be ignored. 

We sample a value of x for each of the 3 A auxiliary field 
variables. Once these values are known, the propagation 
reduces to a rotation in the spin space, and, therefore, 
to multiplying the current spinor value for each particle 
by the set of matrices given by the transformation above. 
For a given eigenvalue A„ < in Ea. lj^UI) the spin states 
of particle k, \ri' k ) = a' k \ f ) + b' k \ \) will be rotated to the 
new one \r]k) having the following components 



a k = a4(cosh(a„) + sinh(a„)^(fc)) 

+ b' k sm^a n )(r n (k)-ir n {k)), 

b k = 6' fc (cosh(a„) - sinh(a„)^(fc)) 

+ a' fe sinhK)(C(fc) +l ^(fc)), (23) 



a n = At\K\^m(k)) 2 + WUW + (^W) 2 , (24) 

and x n is the sampled Hubbard-Stratonovich value. 
For positive values of X n , one has a similar set of equa- 
tions, in which sinh(a„) is substituted with zsin(— a n ). 

Finally it is worth mentioning here that importance 
sampling has been used for the integral in the HS vari- 
ables. The value of the overlap of the walker with the 
trial function will not be peaked around x n = 0, but will 
be shifted. Rather than sampling from the gaussian we 
preferentially sample values where we predict the trial 
function will be larger. One way is to shift the sampled 
gaussian values with a drift term analogous to the drift 
term in diffusion Monte Carlo by replacing the a opera- 
tors by their expectation value at the current i?, S value 
and taking the real part. That is we write 



1 



'2-K J-oo 
1 



dx n e 2 e 



x n ^-\ n AtO n 



/2tt 



dxne- L ^ J ^e^^ Mo -e 



3j ^ — JRjO 
(On) = " 



V-KAt (On) 

y T \O n \R, S) 



(V T \R,S) ' 



(25) 



and sample the shifted gaussian, the last correction term 
is included in the weight. With this real shift and the 
compensating weight, only the efficiency of the algorithm 
is changed. We have tried other schemes using a dis- 
cretized gaussian integration with altered probabilities 
and compensating weights with very little difference in 
the overall efficiency. In Ref. |36( a complex drift rather 
than the real drift in Eq. (|25|l has been used. Unlike our 
real drift above, this can change how the path constraint 
is applied. 



A. Three— body potential 

For a neutron system the spin-dependent part of Ur- 
bana IX potential, given in eqs. Q and JSJl reduces to 
a sum of terms containing only two-body spin operators 
but with a form and strength that depends on the posi- 
tions of three particles. As will be seen below, for a fixed 
position of the particles, the inclusion of three-body po- 
tentials of the Urbana IX type in the Hamiltonian does 
not add any additional complications. It simply changes 
the strength of the coefficients of the terms in the po- 
tential and can be trivially incorporated in the AFDMC 
calculations. 

The anticommutator in Eq.JSJ can be written as 



r, UU 11 i/ 
2 X jki a j a k 



(26) 
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where 



and 

j/ji = Y(m-K,c 3 ,rji) - T(nv,C3,rji), 

= 3T(77v,c 3 ,r j7 )f^ . (28) 

The spin-dependent part of the three-body interac- 
tion V^ SD can then be easily incorporated in the matrix 
Aj t a,k,0 of Ea.lfT5]l. by the following substitution 



.4 



Aj t<X ;k,0 + 2 E B 27r . 



J jkl 



(29) 



Similarly the new terms in the Illinois potentials can 
be included into this matrix. 



B. The Spin-Orbit Propagator 



At this point, Pls is expanded by using the second 
form of this propagator given in Eq. (|31[) keeping both 
linear and quadratic terms in Ar. The integral in R' can 
be done by taking into account that i) the gaussian in- 
tegrates to one if there are no powers of Ar; ii) terms 
containing only one power of a Af integrate to zero; iii) 
quadratic terms containing powers of different compo- 
nents of AR' , integrate to zero and iv) terms like (Ax'j) 2 

integrate to Ath 2 /m. 

We first consider the part coming from the linear terms 
in Ar in both the wave function and Pls- These terms, 
after integration give 



At E 



4i 



[{a 3 +a k )xf ]k ]-W^(R). 



(33) 



The expression above can be further simplified by in- 
terchanging the dummy indices j and k 



A first order approximation |2l| to the spin-orbit con- 
tribution to the propagator can be obtained by operating 
the derivative appearing in the Lj k ■ Sj k operator on the 
free propagator Go 

(Vi-Vfc) G (R,R') = 

717 

-75— (A^ - Ar k )G (R, R!) , (30) 
n At 

and substituting this expression back into the propa- 
gator. As a result, the spin-orbit part Pls of the prop- 
agator is factored out and is finally written as 



P LS = exp E 



mv LS {r jk ) . 
—^!—[r jk x {Ar) jk \-a 3 



exp J2 



r— x r jk) ■ Ar j I (3 1 ) 



K J^k 



where (Ar)j k = Afj — Af k and T,j k = Bj + a k . 

However a careful analysis of the above expressions 
show that they include some spurious contributions lin- 
ear in At. In order to see this the wave function is ex- 
panded, as usual, in the integral form of the imaginary 
time Schrodinger equation keeping only linear terms 



9(R) 



At 



3 



at(R) 



dR'G Q (R, R')P LS MR) 
- ^Ar p -V p *(i?)] + ... 



(32) 



AtJ2v L s(r 3k )[L-S] jk y(R) , 

3<k 



(34) 



which is the spin-orbit contribution to the Hamiltonian 
multiplied by — At. 

However the Pls propagator includes other terms 
which are of the same order in At. They come from 
the quadratic Ar terms of the expansion of Pls 



At(V 2 + V z ) = At— E E E VLs(r jk )v L s(i 
j k ^j pjtj 

(Ej fc x f jk ) ■ (Sj P x f jp ) 

At ^2 E E E v Ls{r 3k )v L s{r 3P ) 

j k ^j p=j Lj 



{fjk 



■ rj p Y,j k ■ Sj p Tijp ■ rj k T,j k ■ rj p 



(35) 



The terms with k = p give rise to a two-body addi- 
tional effective potential V 2 add = — V 2 , 



yadd 



E" 

3<k 



mr % v ls( r jk) 



8h 2 



2 + 3j ■ B k 



crj ■ r jk a k ■ r jk \ 



(36) 



The terms with k ^ p lead to a three-body additional 
effective potential V£ dd = — V3, given by 



yadd 



E E 

j<k<p cyclic 



mrj k rjpV L s{rjk)vLs{rjp) 
16h 2 



{fjk ■ fj p [2 + a k ■ & j + a v ■ Oj + a k ■ B p ] 

3 ' T'jk&k ' T jp @p ' T'jk&j ' 1 jp 
&p ' r jk&k ' T jp } ■ (37) 
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Therefore in the actual propagation it is necessary to 
include explicitly these terms with opposite sign if one is 
using Pls as given by Ea. lpTTJl . 

An alternative method that we have also used comes 
from realizing that the counter terms are produced by 
the next order term in the series expansion of the expo- 
nential. These terms either average to zero, or are higher 
order in the time step or give incorrect contributions. 
Subtracting them gives the propagator, 

mv LS (r jk ) 

eX P I 2^ a.*2 [ r J fc X ( Ar )^] ' °3 



1 



exp 



\ - mv LS (r jk ) . 



(38) 



with the second exponential giving the required counter 
terms to include. The two forms arc equivalent to first 
order in At. 



C. Trial wave function 

In our calculations we use the simple trial function 
given by a Slater determinant of one-body space-spin or- 
bitals multiplied by a central Jastrow correlation, 



j<k 



A 



(39) 



where A is the antisymmetrizer of A particles. The 
overlap of a walker with this wave function is the deter- 
minant of the space-spin orbitals, evaluated at the walker 
position and spinor for each particle (for nuclear matter 
the spinor also includes the isospin), and multiplied by a 
central Jastrow product. 

For unpolarized neutron matter in a box of side L, the 
orbitals are plane waves that fit in the box times up and 
down spinors. The usual closed shells are 2, 14, 38, 54, 
66, 114, ... for neutrons and 4, 28, 76, ... for nucleons. 

The Jastrow correlation function f(r) has been taken 
as the first component of the FHNC/SOC correlation op- 
erator Fij , which minimizes the energy per particle of ei- 
ther neutron or nuclear matter at the desired density |l6j 
(see also Section IV}) . 

As noted in section ITTT1 a trial function with spin ex- 
change and tensor correlations requires exponentially in- 
creasing computational work as the number of particles 
increases. The advantages of our trial function is that it 
is totally antisymmetric and for A particles requires order 
A 3 operations to evaluate. However, it does not contain 
any amplitude generated by the tensor force where spins 
are flipped with a compensating orbital angular momen- 
tum. It is left to the AFDMC method to generate these 
missing components. 



Other forms of a trial wave function can be used. For 
example including a linear combination of Slater deter- 
minants is possible as is modifying the orbitals to in- 
clude spin correlations of backflow form |37|- Both of 
these avoid the exponential computational complexity, 
but may not capture the essential physics of the tensor 
force 381. 



D. Path Constraint 

As in standard fermion diffusion Monte Carlo, the 
AFDMC method has a fermion sign problem. The over- 
lap of our walkers with the trial function will be complex 
in general so the usual fermion sign problem becomes a 
phase problem. 

To deal with this problem, we constrain the path of 
the walkers to regions where the real part of the overlap 
with our trial function is positive. We have also tried 
constraining the phase to that of the trial function as 
in the fixed phase approximation |39j |. Both give about 
the same results, within error bars, and we report values 
where the real part is positive. For spin independent time 
reversal invariant potentials both reduce to the fixed- 
node approximation. It is straightforward to show that 
if the sign of the real part is that of the correct ground 
state, we will get the correct answer and small deviations 
give second order corrections to the energy. We have not 
been able to prove that this constraint always gives an 
upper bound to the ground state energy although it ap- 
pears to do so for the calculations we have done to date. 
It seems likely that there is not an upper bound theorem 
for the mixed estimate of the energy. If forward walking 
or a path integral ground state technique 0, EH is used, 
the method simply produces a better trial function and 
the energy must be an upper bound. 

In the fixed node method 0] the nodal structure of the 
trial function is determined by the Slater Determinant. 
Similary, our path constraint is fully determined by the 
space spin Slater Determinant of Eq. I|39|l . The Jastrow 
function therefore affects only the variance and not our 
final results. 



E. Tail corrections 

Monte Carlo calculations are generally performed 
within the sphere of radius L/2, where L is the length 
of the box side. Usually tail corrections are estimated 
by integrating out the spin-independent part of the two- 
body potential from L/2 up to infinity. We have made 
our calculations within the full simulation box, and, in 
order to include also the contribution from the neighbor 
cells, we have tabulated the Jastrow factor /(r) and the 
components v p (r) of the two-body potential in the fol- 
lowing form 
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F(x,y,z) = JJ f(\(x + mL x )x 

mno 

+ (y + nL y )y + (z + oL z )z\) 
V p (x,y,z) = v p (\(x + mL x )x 

mno 

+ (y + nL y )y+{z + oL z )z\). (40) 

For the calculations shown, we found it adequate to in- 
clude only the 26 additional neighbor cells corresponding 
to m, n, and o taking the values —1, 0, and 1. 

Our results are therefore already tail corrected. We 
found that the standard way of treating tail correc- 
tions leads to results very close to ours, except when we 
consider model interactions which contain tensor forces, 
which are relatively long range forces. 

The three body potential is not treated as the two body 
one. Here we have estimated the tail corrections to the 
three body potential from the PBFHC variational results 
described in|V| This analysis shows that such corrections 
are already very small for systems with 66 nucleons. 



F. The Algorithm 

Finally, in this subsection we give the schematic struc- 
ture of the AFDMC algorithm. 

1. Sample \R, S) initial walkers from \{W T \R, S)\ 2 us- 
ing Metropolis Monte Carlo. 

2. Propagate the spatial degrees of freedom in the 
usual diffusion Monte Carlo way with a drifted 
gaussian for half a time step. 

3. For each walker, diagonalize the potential matrix 
(two- and three-body terms). 

4. Loop over the eigenvectors, sampling the corre- 
sponding Hubbard-Stratonovich variable and up- 
date the spinors for half a time step. Introduce 
approximate importance sampling of the Hubbard- 
Stratonovich variables, as discussed in the previous 
subsection. 

5. Propagate the spin-orbit, using importance sam- 
pling. 

6. Repeat 2, 3, and 4 in the opposite order to produce 
a reversible propagator to lower the time step error. 

7. Combine all weight factors and evaluate the new 
value of (*S?t\R, S). If the real part is less than in- 
clude the walker in the evaluation of the mixed and 
the growth energies, and then enforce constrained 
path by dropping the walker. In general, the im- 
portance sampling makes the number of dropped 
walkers small. 



8. Evaluate the averages of (^>t\R, S), and 
(^>t\H\R, S) to calculate the mixed energy. 

9. Repeat as necessary. 



V. FHNC AND PBFHNC CALCULATIONS 

In this Section we present the method that we have 
used to estimate the finite size effects in AFDMC sim- 
ulations. Such a method is made necessary by the fact 
that simulations with more than 100 nucleons are com- 
putationally very demanding. A many-body theory, like 
FHNC, based on integral equation techniques, in which 
the number of particles in the simulation box has no prac- 
tical limitation seems to be the best candidate to do this. 

FHNC theory was originally developed [4^ to treat 
fcrmionic systems in the thermodynamic limit. However, 
FHNC has been recently reformulated to deal with a fi- 
nite number of fermions in a periodic box, as those used 
for the Monte Carlo calculations |2{| . Such a theory, de- 
noted as Periodic Box FHNC (PBFHNC), is based upon 
the fundamental property of the FHNC cluster expan- 
sion to be valid at all 1/A order (52, an( ^ ^ nas been 
developed for Jastrow-correlated wave functions. In the 
cases of a nucleonic system interacting via a central po- 
tential it has been shown that finite size effects are (i) 
not limited to the kinetic energy expectation value, and 
rather accurately estimated by PBFHNC calculations 

However, realistic correlations F(ij) are spin- 
dependent and have an operatorial structure similar to 
that of the two-body potential, as in Eq.J2Jl (where the 
component p = 1 corresponds to the Jastrow correla- 
tion). Therefore the PBFHNC developed in Ref. 
cannot be used as such, but has to be generalized to 
treat spin-dependent correlations. The main problem is 
that the spin operators involved do not commute, namely 
[F(ij), F(ik)} ^ 0. This feature makes a full FHNC sum- 
mation impossible and one has to resort to reasonable 
approximations for the spin-dependent correlations. 

Such approximations are characterized by the fact 
that, whereas the cluster diagrams containing scalar cor- 
relations only are summed up with FHNC technique, 
only a limited set of cluster diagrams containing spin- 
dependent correlations are included in the calculation. 
The most tested and used approximation is the so called 
FHNC/SOC, described in Ref. 0. In our calculations 
we have used the version adopted in [l6l | in order to 
compute the different correlation functions F(i,j) at the 
various densities considered. We have considered only 
the three variational parameters corresponding to heal- 
ing distance, d c , of central (p = 1 — 4) and spin-orbit 
correlations (p = 7,8), the healing distance, d t of tensor 
correlations (p = 5,6), and the quencher, a s , of the spin- 
isospin dependent correlation. The other variational pa- 
rameters, like the spin-independent potential quencher 
and the correlation quenchers have been kept fixed at 
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TABLE I: Variational parameters used in our FHNC/SOC 
and PBFHNC calculations for the AUG' and AU8' poten- 
tials, ro = (3/(47rp)) 1/ ^ 3 is the average distance between the 
neutrons, ro, d c and dt are given in fm. The reference density 
po = 0.16 fm~ 3 is the equilibrium density of nuclear matter. 



p/po 


ro 


d c (6) 


dt(6) 


a s (6) 


d c (8) 


*(8) 


a 8 (8) 


0.75 


1.258 


1.761 


4.695 


0.9 


2.264 


4.528 


0.8 


1.00 


1.143 


1.714 


4.571 


0.9 


2.285 


4.571 


0.8 


1.25 


1.061 


1.485 


4.752 


0.9 


2.228 


3.960 


0.8 


2.0 


0.907 


1.723 


4.595 


0.8 


2.267 


4.535 


0.7 


2.5 


0.842 


1.768 


4.715 


0.8 


2.189 


5.004 


0.7 



TABLE II: FHNC/SOC energy per particle of neutron mat- 
ter for the AU6' interaction, at various densities. Tf is the 
Fermi kinetic energy, and (T) is the kinetic energy expecta- 
tion value, corresponding to the average of the JF and PB 
kinetic energies. (Vjs) and (V3) are the expectation values of 
the two-body and three-body potentials respectively. AE2 is 
the second order perturbative correction 45]. AE e i em is the 
contribution from the lowest order elementary diagram (see 
text). All the quantities, except p/po, are expressed in MeV. 



P/po 


T F 


(T) 


(V2) 


(Vs) 


-Efhnc 


AE 2 




0.75 


28.969 


35.33 


-22.67 


2.58 


15.2 


-0.9 


0.6 


1.00 


35.094 


43.82 


-28.58 


5.17 


20.4 


-0.9 


0.9 


1.25 


40.722 


52.27 


-34.11 


8.53 


26.7 


-1.5 


1.2 


2.0 


55.708 


74.40 


-46.93 


27.29 


54.8 


-4.4 


2.8 


2.5 


64.643 


88.85 


-53.36 


44.72 


80.2 


-6.1 


3.8 



unity. The optimal values of such variational parame- 
ters for pure neutron matter are shown in Tabled They 
have been obtained by minimizing the average energy 
E av = \{Ejf + Epb), where the two energy expecta- 
tion values Ejp and Epb refer to the Jackson-Feenberg 
and Pandharipande-Bethe kinetic energy expressions re- 
spectively 0. The usual constraint |gjJ E ~f FB| < .005 
has been imposed in order to limit the range of vari- 
ability of the free parameters in a region of reliability of 
the FHNC/SOC approximation. We have verified that in 
such region the normalization condition is fulfilled within 
a few percent. 

The variational energies for the case of the AU6' inter- 
action are reported on Table [H] The Table also reports 
the second order CBF perturbative corrections AE 2 
and the contribution from the lowest order elementary 
diagram AE e i em , as discussed in Ref. . The not neg- 
ligible value of AE e i em indicates that the effect from ele- 
mentary diagrams is larger than has been assumed in all 
the past FHNC/SOC calculations of the nuclear matter 
equation of state H H|. In recent FHNC/SOC cal- 
culations of the equation of symmetric nuclear matter 
and pure neutron matter |18l l4q extra cluster diagrams 
with respect to the approximation used here have been 
included. Differences between the various FHNC/SOC 
calculations are within the predictive accuracy of the ap- 
proximation. 

In Tabic IIIII we compare the results of two different 



TABLE III: Comparison of the FHNC/SOC results for the 
AU8' interaction, obtained with correlation operator of the 
type fa or of the type /8. In the first case the contribution of 
the spin-orbit potential is calculated perturbatively from the 
AUG' Hamiltonian. For comparison, in the third column the 
results for the AU 6' interaction are also reported. In all cases 
the contribution from elementary diagrams has been added. 



p/po 


T F 


AUG' 


AU8'(f 6 ) 


AU8'{h) 


0.75 


28.969 


15.8 


16.1 


13.3 


1.00 


35.094 


21.3 


21.8 


17.6 


1.25 


40.722 


27.9 


28.8 


23.0 


2.0 


55.708 


57.6 


59.0 


47.5 


2.5 


64.643 


84.0 


86.2 


71.7 


TABLE IV: Comparison of the energy E2 at the second order 


of the FHNC cluster expansion 


with the full FHNC energy, 


-BpBFHNC- 


The calculation has been performed for the v' 6 


model interaction at p 


= 0.16 fm and Jastrow correlation 


factor. 










N 


T F 




E 2 


-EpBFHNC 


14 


35.600 




19.36 


17.60 


38 


33.703 




17.51 


15.91 


66 


34.917 




19.11 


17.63 


114 


35.646 




20.09 


18.71 


1030 


35.139 




19.46 


18.04 



FHNC/SOC calculations of the equation of state of neu- 
tron matter, carried out for the AU8' potential. In the 
first one (AU8' (fa)) the spin-orbit correlation is set equal 
to zero, whereas, in the second one (AU8'(f 8 )), is in- 
cluded. One can see that the introduction of the spin- 
orbit correlation leads to a large lowering in the energy. 
As it will be shown, we do not find such a lowering when 
the spin-orbit interaction is included in the AFDMC sim- 
ulations. In the FHNC/SOC approximation the cluster 
contributions from spin-orbit correlations are correctly 
included only at the lowest order level. The many- 
body cluster contributions are essentially neglected. The 
large and attractive spin-orbit contribution found in the 
AU8'(fs) calculation may be due to this inaccuracy. On 
the other hand it might be possible that nodal surface 
induced by the spin-orbit part of the interaction is not 
accurately described by our trial function. 

In order to compute the finite size effects in a realistic 
way one should first generalize the PBFHNC theory to in- 
clude SOC diagrams like in FHNC/SOC approximation. 
Work in this direction is in progress [47] . In this paper 
we limit ourselves to including only the two-body cluster 
diagrams for the two-body potential and the kinetic en- 
ergy and the leading three-body cluster diagrams for the 
three-body potential H3 in the PBFHNC scheme. Such 
leading terms correspond to include up to two correla- 
tion operators in the three-body cluster diagrams. We 
will show that this approximation, hereafter denoted as 
PBFHNC/L, can already be used to roughly estimate the 
finite size effects. 



10 



TABLE V: PBFHNC/L results for the AU6' interaction at 
density p — 0.16fm -3 . The Fermi kinetic energy Tf, the ex- 
pectation values of the kinetic energy (T), the two-body po- 
tential (V)2 and the three-body potential (V)3 are displayed 
together with the energy per particle E in MeV units. 



N 


T F 


<T> 


(V)2 


(V) 3 


E 


14 


35.600 


44.47 


-29.41 


4.31 


19.37 


38 


33.703 


42.41 


-29.43 


4.70 


17.68 


66 


34.917 


43.64 


-29.07 


4.82 


19.39 


114 


35.646 


44.40 


-28.87 


4.87 


20.40 


1030 


35.139 


43.88 


-28.95 


4.86 


19.79 


TABLE VI: As 


in Table IVI at density p 


= 0.32 far 


-3 


N 


T F 


(T) 


(V)2 


(Vh 


(E) 


14 


56.512 


74.33 


-48.04 


17.18 


43.47 


38 


53.500 


71.64 


-50.25 


19.36 


40.75 


66 


55.428 


73.41 


-49.51 


20.30 


44.20 


114 


56.584 


74.56 


-48.94 


20.78 


46.40 


1030 


55.779 


73.75 


-49.08 


20.84 


45.51 



The performance of the two-body cluster approxima- 
tion to account for finite size effects is studied in Table llVl 
There, for a purely central potential without three-body 
force, PBFHNC/L and PBFHNC energies are compared 
at p — 0.16 fm -3 for the range of particle numbers used 
in our quantum Monte Carlo simulations. 

Tables M and IVTl give the PBFHNC/L results for the 
AU6' interaction at two different densities for a number 
of neutron systems. Notice that the energy differences 
between the cases with 66 and 114 neutrons are very close 
to those obtained in the AFDMC simulations, given in 
Tables lVIIII and lrXl Systems with 14 and 38 neutrons are 
too small to be included in the finite size effects analysis. 

VI. RESULTS 
A. AFDMC results for neutron matter 

Extensive neutron matter calculations have been car- 
ried out for the AU6' and AU8' interactions by consid- 
ering 14, 38, 66 and 114 neutrons in a periodic box at 
various densities ranging from 0.75po up to 2.5po- 

In figure ^ we show a typical behavior of the mixed 
and growth energy as a function the time step for 14 
neutrons in a periodic box at p = 0.32 fm -1 interacting 
via AU8'. At At = 5 x 10~ 5 far 1 we have found that 
the statistical error are smaller than the extrapolation 
ones irrespective of the density and number of particles. 
All the calculations reported here have been obtained by 
using that value for the time step. 

The 14 neutrons system is interesting because it is 
small enough to be studied by using other many body 
methods which become inefficient for larger systems. In 
order to provide a full set of results for this system in 



50 I 1 1 1 1 1 1 1 — 

mixed 




6c-0o 0.0001 0.00014 0.00018 

At (MeV" 1 ) 

FIG. 1: Mixed and Growth energies versus the time step for 
14 neutrons at p = 0.32 fm -3 with the AU8' interaction. 



TABLE VII: AFDMC energies per particle in MeV of 14 
neutrons in a periodic box for interaction models at various 
densities. Error bars for the last digit are shown in parenthe- 
ses. 



P(fm 3 ) Ve ^8 

0.12 12.41(4) 12.32(5) 

0.16 15.12(4) 14.98(6) 

0.20 17.86(5) 17.65(7) 

0.32 27.84(6) 27.3(1) 

0.40 36.0(1) 35.3(1) 



table IVIll we report the energies at several densities cal- 
culated with the v' e and v' s interactions. 

Diffusion Monte Carlo calculations using a pair- 
product wave function for 14 neutron systems have just 
been reported [28|. They however set the potential dis- 
continuously to zero at distances greater than L/2, while 
we use either the nearest image convention or a lattice 
sum giving a continuous potential. We expect better ex- 
trapolation to large system sizes with the continuous po- 
tential as well as smaller time step errors. The time step 
errors will affect our AFDMC calculations more because 
we currently use the primitive approximation rather than 
building the Green's function from a product of exact 
two-body Green's functions. In principle we could use 
the Hubbard-Stratonovich breakup for the pair-product 
Green's function. In any case, we have carried out a cal- 
culation at p = 0.16 fm~3 using the same discontinuous 
potential and obtained 20.64(2) MeV and 20.32(6) MeV 
for the Vq and v' s potentials respectively compared with 
their values of 19.91(11) and 17.00(27). The larger differ- 
ence when the spin-orbit term is included in the Hamil- 
tonian may be due to the different trial wave functions 
used. 

In Tables IVTTll and llXl we report the results obtained 
with the AFDMC method of this work for neutron matter 
at the different densities considered for various system 
sizes. The extrapolation to infinite number of particles is 
carried out by using the PBFHNC/L results for a given 
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0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 

p (far 3 ) 

FIG. 2: AFDMC energy per particle for neutron matter with 
obtained from simulations with 66 neutrons and the v' s poten- 
tial. The variational FHNC/SOC results obtained with cor- 
relation functions of type Fe and Fg are and the Brueckner- 
Hartree-Fock (BHF) results of also plotted and the The 
statistical error in the AFDMC results are smaller than the 
symbols. 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 

p (far 3 ) 

FIG. 3: Extrapolated AFDMC Equation of state of pure 
neutron matter with the AU8' potential (solid line). The 
variational results of Refs. |4(J (APR, dotted-dashed line) 
and fist (MPR, dashed line) corresponding to the Argonne 
vis -two body plus Urbana IX -three body potential are also 
plotted. The lines are for guiding the eyes. The statistical 
errors of the AFDMC estimates are smaller than the symbols. 



number of neutrons and for the infinite system. 

The spin-orbit contribution is rather small at all of 
the densities considered. This contrasts with previous 
FHNC/SOC calculations. In Fig. we plot the AFDMC 
results together with the variational FHNC/SOC results 
for the v' s interaction obtained by using correlation oper- 
ators of the Fq and F$ forms, and the Brucckncr-Hartree- 
Fock (BHF) results for the vis potential [48j ■ One can see 
that SOC(F6) and SOC(F8) in the figure give quite dif- 
ferent equations of state, particularly at high density. We 
have tried transient estimates and AFDMC simulations 
with orbitals of spin-backflow type for the 14 neutrons 
system finding not more than roughly 5% lowering of the 
energy |3S|. 

The three-body potential gives a large contribution to 
the energy per particle at high densities. Therefore the 
search for a realistic three-body potential is a very funda- 
mental problem for the study of dense and cold hadronic 
matter. A considerable amount of work has been done 
to find, in a semi-phcnomcnological way, three-body po- 
tentials to describe ordinary matter. However, whether 
such potentials are also valid in the high density regime 
is still an open and debated issue. In table Ixl we report 
AFDMC results performed with the two body v' 6 inter- 
action and five different three-body potentials including 
the Urbana IX [T3, E3| (UEK) and the recent Illinois 1 
through 4 0] three-body interactions. One can see that 
already at twice the nuclear matter density, the energy 
contributions from the three-body potentials are large 
and very different from each other, in spite of the fact 
that all of them provide a satisfactory fit to the ground 
state and the low energy spectrum of nuclei with A < 8. 

In Fig. |3 we show the AFDMC equation of state for 
pure neutron matter with the AU8' interaction corre- 
sponding to the extrapolated values for infinite matter. 
We compare with the variational results of Akmal et al. 



[iq and the more recent ones of Ref. [Tg] both of them 
obtained with the Argonne wis 2- and Urbana IX 3- nu- 
cleon interactions. One can see that there is a surprising 
good agreement between our AFDMC results and the 
latest variational calculation of Morales, Jr. et al. [l8| . 
The compressibility /C, given by 



1 _ p3 d 2 E (p) 



dp 2 



2p 



2 dE (p) 



dp 



(41) 



can be estimated from the equation of state by tak- 
ing Eq — E/A. For a Fermi gas the compressibility is 
K. F = 9Tr 2 m/(k 5 f h 2 ). The AFDMC results for JC/K F 
obtained from the extrapolated energies with AU8' are 
shown in Fig. 0] They are compared with the compress- 
ibility calculated from the variational energies of Ref. 

Mm. 



B. AFDMC results for nuclear matter 

The AFDMC can deal with N ^ Z systems, and we 
have applied it to compute the asymmetry coefficient of 
the mass formula for the semirealistic two-body potential 
MS3 which is spin-isospin dependent but has no tensor 
force |49|,|5(j The resulting values of E/A at po for sym- 
metrical nuclear matter are given in TableEU where they 
are also compared with the FHNC/SOC and PBFHNC 
results. The finite size correction is estimated from the 
corresponding PBFHNC results. 

In Fig. we plot the AFDMC energy per particle 
as a function of the asymmetry parameter, a — (N — 
Z)/(N + Z), of nuclear matter. The FHNC/SOC curve 
corresponds to a quadratic fit of nuclear matter (a = 0) 
and pure neutron matter (a = 1). 

FHNC/SOC can only be used to study N = Z or 
N = A matter. The symmetry energy obtained from 
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TABLE VIII: AFDMC energies per particle in MeV for the AU& interaction obtained with systems with 14, 38, 66 and 114 
neutrons at various densities. Error bars for the last digit of the Monte Carlo calculations are shown in parentheses. The last 
column gives the extrapolated values from the PBFHNC/L calculation |47f . 



/o(fm~ 3 ) AFDMC(14) AFDMC(38) AFDMC(66) AFDMC(114) AFDMC(oo) 

0.12 14.96(6) 13.76(9) 14.93(4) 15.62(8) 15.0 

0.16 19.73(5) 18.56(8) 20.07(5) 20.99(9) 20.4 

0.20 25.29(6) 24.4(1) 26.51(6) 27.6(1) 26.9 

0.32 48.27(9) 49.8(1) 53.11(9) 55.3(2) 54.4 

0.40 69.9(1) 74.5(2) 79.4(2) 82.2(2) 81.3 



TABLE IX: AFDMC energies per particle in MeV for the AU8' interaction obtained with systems with 14, 38 and 66 neutrons 
at various densities. Error bars for the last digit of the Monte Carlo calculations are shown in parentheses. The last column 
gives the extrapolated values from the PBFHNC/L calculation [4^ |. 



p(fm~ 3 ) AFDMC(14) AFDMC(38) AFDMC(66) AFDMC(oo) 

0.12 14.80(9) 13.96(5) 15.26(5) 15.5 

0.16 19.76(6) 18.67(6) 20.23(9) 20.6 

0.20 25.23(8) 24.7(1) 27.1(1) 27.6 

0.32 48.4(1) 46.8(2) 54.4(6) 55.6 

0.40 70.3(2) 76.3(2) 81.4(3) 83.5 



TABLE X: AFDMC energies per particle in MeV for the 
t>6+IL potentials calculated with 66 particles. For the case of 
i>e+IL2 interaction, at p = 0.32 fm -3 the energies per par- 
ticle with 38 and 54 neutrons are 12.6(2) and 10.0(3) MeV 
respectively. 



p(fm" 3 ) AU6' 


IL1 


IL2 


IL3 


IL4 


0.16 20.07(5) 


11.2(1) 


11.39(8) 


12.0 (4) 


10.5(2) 


0.32 53.11(9) 


8.0(4) 


11.1(3) 


14.7(3) 


4.7(3) 
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FIG. 4: Compressibility ratio, K./K.F, for neutron matter 
obtained from the extrapolated AFDMC energy per particle 
with the AU8' potential (solid line). The compressibility ob- 
tained from the variational results of Refs. |4q| (APR, dotted- 
dashed line) and [is| (MPR, dashed line) is also plotted. The 
lines are for guiding the eyes. The statistical errors of the 
AFDMC estimates are smaller than the symbols. 



FHNC/SOC is 41.59 MeV. The function E/A(a) pro- 
vided by the AFDMC results is not fully quadratic in a, 
and corresponds to a symmetry energy of ~ 36.4 MeV. 



TABLE XI: Finite size corrections for symmetrical nuclear 
matter PBFHNC results for the MS3 potential at 

p = 0.16 fm" 3 . The PBFHNC calculations have been per- 
formed with a Jastrow correlated wave function, whereas the 
FHNC/SOC result has been obtained with a correlation oper- 
ator of the type F 4 . PBFHNC and FHNC/SOC calculations 
include the basic four-point elementary diagram E4. 



A 


PB-FHNC 


FHNC/SOC 


AFDMC 


28 


-13.6 




-16.17(6) 


76 


-15.6 




-18.08(3) 


2060 


-14.0 




1 


00 


-14.0 


-14.9 


-16.5(1) 



VII. OUTLOOK AND CONCLUSIONS 

We have described a quantum Monte Carlo method 
specially suited to perform calculations on nucleon sys- 
tems with noncentral interactions. It has been applied 
here to calculate the equation of state of pure neutron 
matter with fully realistic interactions by approximating 
it with up to 114 neutron in a simulation box. Finite size 
effects have been estimated by performing 2- plus 3- body 
cluster diagrams calculations based on PBFHNC method 
with spin-dependent correlations. The results obtained 
show an overall agreement with Brueckner-Hartree-Fock 
calculations and with a recent 2- plus 3- body cluster di- 
agrams variational calculation |18|. They also indicate 
that there is a very small contribution coming from the 
spin-orbit component of the two-body interaction while 
the effect from the three-body potential is quite large, 
particularly at high densities. The large differences ob- 
tained for the equation of state for different phenomeno- 
logical three-body potentials point out a three-body po- 
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FIG. 5: AFDMC and FHNC/SOC energy per particle of 
nuclear matter for several values of the asymmetry parameter 
[25| . The lines correspond to polynomial fits of the calculated 
energies. 



We believe that our method should be able to pro- 
duce accurate Monte Carlo calculations of a wide variety 
of nuclear systems. While previous Monte Carlo calcula- 
tions have been severely restricted on the particle number 
by the spin-isospin sum this restriction is lifted by using 
the auxiliary field breakup of the spin-isospin part of the 
Hamiltonian, while using standard diffusion Monte Carlo 
for the spatial degrees of freedom. A pressing need is sim- 
ulating nuclear matter with fully realistic interactions as 
already done for pure neutron matter; calculating the 
properties of light nuclei to compare with exact GFMC 
calculations; and investigating pion condensation. In ad- 
dition, including explicit meson degrees of freedom can 
also be attempted. In the language of this paper, each 
meson field mode corresponds to an auxiliary field |22j |. 



tential problem in the study of dense and cold hadronic 
matter. 

Work in progress is to validate the present results using 
trial wave functions, other than the simple Slater deter- 
minant given in Eq. JSHJ, that also can be calculated 
efficiently j3^| ■ This will allow us to both lower the vari- 
ance of our calculations, as is usual when better guiding 
functions are used in the importance sampling of the ran- 
dom walk, as well as to obtain a better path constraint. 



Acknowledgments 

We wish to thank A. Fabrocini, V. R. Pandharipande, 
A. Polls, S. Pieper, and R. Wiringa for helpful conver- 
sations. S. F. wish to thank the International Centre 
for Theoretical Physics in Trieste for partial support. 
A. S. acknowledges the Spanish Ministerio de Ciencia y 
Tecnologia for partial support under contract BMF2002- 
00200 



[1] V. R. Pandharipande, I. Sick, and P. K. A. deWitt Hu- 
berts, Rev. Mod. Phys. 69, 981 (1997). 

[2] G. G. Raffelt, The stars as laboratories of fundamen- 
tal physics (University of Chicago, Chicago & London, 
1996). 

[3] R. F. Sawyer, Phys. Rev. D 11, 2740 (1975). 
[4] R. F. Sawyer, Phys. Rev. C 40, 865 (1989). 
[5] N. Iwamoto and C. J. Pethick, Phys. Rev. D 25, 313 
(1982). 

[6] S. Reddy, M. Prakash, J. M. Lattimer, and J. A. Pons, 

Phys. Rev C 59, 2888 (1999). 
[7] R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla, Phys. 

Rev C 51, 38 (1995). 
[8] V. G. J. Stoks, R. A. M. Klomp, C. P. F. Terheggen, and 

J. J. de Swart, Phys. Rev C 49, 2950 (1994). 
[9] R. Machleidt, F. Sammarruca, and Y. Song, Phys. Rev 

C 53, R1483 (1996). 
[10] V. G. J. Stoks, R. A. M. Klomp, M. C. M. Rentmeester, 

and J. J. de Swart, Phys. Rev C 48, 792 (1993). 
[11] L. Engvik, M. Hjorth- Jensen, R. Machleidt, H. Miither, 

and A. Polls, Nucl. Phys. A 627, 85 (1997). 
[12] B. S. Pudliner, V. R. Pandharipande, J. Carlson, and 

R. B. Wiringa, Phys. Rev. Lett. 74, 4396 (1995). 
[13] S. C. Pieper, V. R. Pandharipande, R. B. Wiringa, and 

J. Carlson, Phys. Rev. C 64, 14001 (2001). 
[14] S. C. Pieper, K. Varga, and R. B. Wiringa, Phys. Rev. 

C 66, 044310 (2002). 
[15] A. Ramos, W. H. Dickhoff, and A. Polls, Phys. Rev. C 

43, 2239 (1991). 



[16] R. B. Wiringa, V. Fiks, and A. Fabrocini, Phys. Rev. C 

38, 1010 (1988). 
[17] A. Akmal and V. R. Pandharipande, Phys. Rev. C 56, 

2261 (1997). 

[18] J. Morales, Jr., V. R. Pandharipande, and D. G. Raven- 
hall, Phys. Rev. C 66, 054308 (2002). 

[19] K. E. Schmidt and M. H. Kalos, in Monte Carlo Meth- 
ods in Statistical Physics, edited by K. Binder (Springer- 
Verlag, Berlin, 1984), pp. 125-143. 

[20] B. S. Pudliner, V. R. Pandharipande, J. Carlson, S. C. 
Pieper, and R. B. Wiringa, Phys. Rev C 56, 1720 (1997). 

[21] S. C. Pieper, in Microscopic Quantum Many-Body The- 
ories and their applications. Lecture Notes in Physics, 
edited by J. Navarro and A. Polls (Springer- Verlag, 
Berlin, 1998), vol. 510, p. 337. 

[22] K. E. Schmidt and S. Fantoni, Phys. Lett. B 446, 93 
(1999). 

[23] S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. 
Lett. 74, 3652 (1995). 

[24] S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. 
B 55, 7464 (1997). 

[25] S. Fantoni, A. Sarsa, and K. E. Schmidt, in Advances in 
Quantum Many-Body Theory, edited by R. F. Bishop, 
K. A. Gernoth, and N. R. Walet (World Scientific, Sin- 
gapore, 2001), vol. 5, pp. 143-151. 

[26] S. Fantoni, A. Sarsa, and K. E. Schmidt, Prog. Part. 
Nucl. Phys. 44, 63 (2000). 

[27] S. Fantoni, A. Sarsa, and K. E. Schmidt, Phys. Rev. Lett. 
87, 181101 (2001). 



14 



J. Carlson, J. Morales, Jr., V. R. Pandharipande, and 

D. G. Ravenhall (2003), nucl-the/0303041. 

S. Fantoni and K. E. Schmidt, Nucl. Phys. A 690, 456 
(2001). 

R. B. Wiringa and S. C. Pieper, Phys. Rev. Lett. 89, 
182501 (2003). 

J. Carlson, V. R. Pandharipande, and R. B. Wiringa, 
Nucl. Phys. A 401, 59 (1983). 

J. Lomnitz-Alder, V. R. Pandharipande, and R. A. 
Smith, Nucl. Phys. A 361, 399 (1981). 
J. Carlson, Phys. Rev. C 38, 1879 (1988). 

E. Buendi'a, F. J. Galvez, J. Praena, and A. Sarsa, J. 
Phys. G: Nucl. Part. Phys. 26, 1795 (2000). 

J. Carlson and M. H. Kalos, Phys. Rev. C 32, 2105 
(1985). 

S. Zhang and H. Krakauer (2003), cond-mat/0208340. 
S. Fantoni, in Advances in quantum many body theories, 
edited by A. Fabrocini, S. Fantoni, and E. Krotscheck 
(World Scientific, Singapore, 2002), vol. 7, pp. 379-405. 
L. Brualla, S. A. Vitiello, A. Sarsa, S. Fantoni, and K. E. 
Schmidt, in preparation (2003). 

G. Ortiz, D. M. Ceperley, and R. M. Martin, Phys. Rev 
Lett. 71, 2777 (1993). 

S. Baroni and S. Moroni, Phys. Rev. Lett. 82 , 4745 
(1999). 

A. Sarsa, K. E. Schmidt, and W. R. Magro, J. Chem. 
Phys. 113, 1366 (2000). 

S. Fantoni and S. Rosati, Nuovo Cimento A 10, 145 
(1974); Lett. Nuovo Cimento 10, 545 (1974); Nuovo Ci- 
mento A 25, 593 (1975). 

S. Fantoni and A. Fabrocini, in Microscopic Quan- 
tum Many-Body Theories and their applications. Lec- 
ture Notes in Physics, edited by J. Navarro and A. Polls 
(Springer- Verlag, Berlin, 1998), vol. 510, pp. 119-186. 
V. R. Pandharipande and R. B. Wiringa, Rev. Mod. 
Phys. 51, 821 (1979). 

A. Fabrocini, private communication (2002). 

A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, 

Phys. Rev. C 58, 1804 (1998). 

F. Arias de Saavedra, A. Sarsa, S. Fantoni, and K. E. 
Schmidt, in preparation (2003). 

M. Baldo, G. Giansiracusa, U. Lombardo, and H. Q. 
Song, Phys. Lett. B 473, 1 (2000). 



[49] I. R. Afnan and Y. C. Tang, Phys. Rev. 175, 1337 (1968). 

[50] R. Guardiola, in Recent progress in Many-Body Theories. 
Lecture Notes in Physics, edited by J. G. Zabolitzky, 
M. de Llano, M. Fortes, and J. W. Clark (Springer- 
Verlag, Berlin Heidelberg, 1981), vol. 142, pp. 398-406. 



APPENDIX A: POTENTIAL PARAMETERS 



The A p _ m parameters in the expansion of the Argonnc 



potentials, Eq. are shown in Table IXIII 

The F m functions in Eq. are written as follows 



Fi(r) 


= T 2 (fi,c;r) , 




F 2 (r) 


= (1 + a r)W(r) 


j 


F 3 (r) 


= firW(r) , 




Fi{r) 


= (^r) 2 W(r) , 




F 5 (r) 


= aiY(m ,c;r) - 


a,2rW(r) 


F 6 (r) 


= a 3 Y(m c ,c;r) - 




F 7 (r) 


= aiT(m ,c;r) , 




Fs(r) 


= a 3 T(m c ,c;r) , 





(Al) 



where the tensor, Yukawa and Wood-Saxon functions are 
defined as 



T(m, c; r) 

Y(m, c; r) 
W(r) 



mr 



"(I 



(mr) 2 ' 



"(I 



1 + cxp( 5(r- 0.5)) 



(A2) 



The the coefficients ao, ■ ■ ■ , »4 are shown in Table IXTTTl 
and the masses mo , m c and /i and the cut-off parameter 
c are given in Table IxTvl 

The values of the parameters of the three-body Urbana 
IX potential are shown in table IX VI 
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TABLE XII: Argonne v8' two-body potential. Matrix A p>m appearing in Eq.@. 



p 


A , 


A n 


A o 


A . 
Sip, 4 


A r 


A ^ 


A » 


A „ 


1 

1 


*7 7/11 

- 1 .o22ol (41 


in i p. onno An An 
2616.39024949 


U 


14/ . (9390526 


U 


U 


U 


u 


2 


-0.12318501 


84.20118403 





-61.22868919 














3 


0.48726001 


-82.48240972 





49.26463509 














4 


0.65399916 


-107.98800762 





-20.40956306 


1/3 


2/3 








5 


0.94963459 


-2.91931242 


-424.28015518 


-398.23289299 














6 


-0.17865545 


-0.97310414 


234.18526077 


-256.12175941 








1/3 


2/3 


7 


-0.71193373 


-373.43774331 





653.08534247 














8 


-0.28568125 


-201.79028547 





354.25604242 















TABLE XIII: Argonne v8' two-body potential. Values of the 
strength parameters in eqs. 1AH and lA2ll 

ao (fm -1 ) ai a^fm -1 ) as 04 (fm~ ) 

0.37929090 3.15588245 10.48427302 3.48918764 11.21004425 



TABLE XIV: Argonne v8' two-body potential. Values of the 
masses and cut-off parameters in eqs. 1AH and lA2ll 

mo (fm _1 ) m c (fm _1 ) /^(fm -1 ) c(fm -2 ) 

0.68401113 0.70729025 0.69953054 2~X 



TABLE XV: Urbana IX three-body potential. Values of the 
parameters appearing in Eq.(JHJ. 



(MeV) 


U (MeV) 


m w (fm 1 ) 


c 3 (far 2 ) 


-0.0586 


0.0048 


0.69953054 


2.1 



